Competing Risks

Rates, Probabilities, and Embedding

Learning Objectives

Pitfalls of traditional transition probability conversion formulas.
Discuss compound transitions and how to account for them.
Demonstrate how to properly embed a transition probability matrix into a timestep.

Code

source("./manifest.r")

# # Better accuracy that "life-table" aka trapezoidal method
# alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8), 49, 43, 59, 17) / 48
# alt_simp      <- function(x,h) h*sum(alt_simp_coef(length(x)) * x)
# 
# r_H_CVD <- 0.15
# r_HD <- .006
# hr_S <- 10
# r_SD <- hr_S*r_HD
# cycle <- 1
# 
# # 1-exp(-r_H_CVD*cycle)
# # 1-exp(-r_HD*cycle)
# 
# get_P <- function(r_H_CVD,r_HD,r_SD,cycle) {
#   P_st <- 
#   matrix(
#   c(exp(-r_H_CVD*cycle) + exp(-r_HD*cycle) -1 ,  (1-exp(-r_H_CVD*cycle)) ,  1-exp(-r_HD*cycle), 
#   0, exp(-r_SD*cycle), 1-exp(-r_SD*cycle),
#   0,0,1),
#   nrow = 3, 
#   ncol = 3,
#   byrow=TRUE,
#   dimnames = list(c("Healthy","Sick","Dead"), c("Healthy","Sick","Dead"))
#   ) %>% 
#     data.frame() %>% 
#     rownames_to_column(var = "from") %>% 
#     mutate(name = "P_st")
# 
#   P_eQ <- matrix(
#     c(-(r_H_CVD*cycle+r_HD*cycle),r_H_CVD*cycle,r_HD*cycle,
#       0,-r_SD*cycle,r_SD*cycle,
#       0,0,0),
#       nrow=3,
#       ncol=3,
#       byrow=TRUE,
#     dimnames = list(c("Healthy","Sick","Dead"), c("Healthy","Sick","Dead"))
#   ) %>% expm() %>% as.matrix() %>% 
#     data.frame() %>% 
#     rownames_to_column(var = "from") %>% 
#     mutate(name = "P_eQ")
# 
# 
#   P_cr <- 
#     matrix(
#     c(exp(-(r_H_CVD+r_HD)*cycle) , (r_H_CVD/(r_H_CVD+r_HD)) * (1-exp(-(r_H_CVD+r_HD)*cycle))   ,(r_HD/(r_H_CVD+r_HD)) * (1-exp(-(r_H_CVD+r_HD)*cycle))   , 
#     0, exp(-r_SD*cycle), 1-exp(-r_SD*cycle),
#     0,0,1),
#     nrow = 3, 
#     ncol = 3,
#     byrow=TRUE,
#     dimnames = list(c("Healthy","Sick","Dead"), c("Healthy","Sick","Dead"))
#   ) %>% 
#     data.frame() %>% 
#     rownames_to_column(var = "from") %>% mutate(name = "P_cr")
#   
#   
#   
#   out <- bind_rows( P_st, P_cr , P_eQ)
#   return(out)
# }

Part 1: A World Without Death

We will first construct a very simple CVD model with only two states: Healthy and CVD.
Transitions are goverened by a single rate from Healthy CVD

We will first construct a very simple CVD model with only two states: Healthy and CVD.
Transitions are goverened by a single rate from Healthy CVD

Parameterize the model

params = 
  list(
    t_names = c("natural_history"),   # Strategy names. 
    n_treatments = 1,                 # Number of treatments
    s_names  = c("Healthy", "CVD"),   # State names
    n_states = 2,                     # Number of states
    n_cohort = 1,                     # Cohort size
    n_cycles = 100,                   # Number of cycles in model.  
    cycle = 1,                        # Cycle length
    initial_age = 55,                 # Cohort starting age
    r_H_CVD = 0.15                    # Rate of healthy -> CVD
  )

Transition Probability Matrix

We’ll next construct a transition probability matrix using the standard rate-to-probability conversion formula:

\[ p = 1 - \exp(-rt) \] where \(r\) is the rate, \(t\) is the time-step (e.g., 1 for annual, 1/12 for monthly, etc.).

Transition Probability Matrix

params$mP <- 
  with(params,{
    tp_H_CVD = 1 - exp(-r_H_CVD * cycle) 
    array(data = c(1 - tp_H_CVD, 0, 
                 tp_H_CVD,1),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
  })

Transition Probability Matrix

params$mP <- 
  with(params,{
    tp_H_CVD = 1 - exp(-r_H_CVD * cycle) 
    array(data = c(1 - tp_H_CVD, 0, 
                 tp_H_CVD,1),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
  })

Next, place these probabilities within the transition probability matrix array.

Outer array dimensions goverened by the number of strategies under consideration.

Transition Probability Matrix

params$mP <- 
  with(params,{
    tp_H_CVD = 1 - exp(-r_H_CVD * cycle) 
    array(data = c(1 - tp_H_CVD, 0, 
                 tp_H_CVD,1),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
  })

Remaining code simply names the various dimensions of the array.

Transition Probability Matrix

$natural_history
         to
from      Healthy     CVD
  Healthy 0.86071 0.13929
  CVD     0.00000 1.00000

Construct the Markov Trace

sim_cohort <- function(params) {
  tr_ <- t(c("Healthy" = params$n_cohort, "CVD" = 0)) 
  trace <- 
    0:params$n_cycles %>% map_df(~({ 
        apply(params$mP,3, simplify=FALSE, function(tp) {  
          (tr_ %*% (tp %^% .x)) %>% data.frame()
        })
    })) %>% 
    as.matrix() 
}

sim_cohort(params) is a function that constructs a Markov trace from the supplied parameters.

Construct the Markov Trace

sim_cohort <- function(params) {
  tr_ <- t(c("Healthy" = params$n_cohort, "CVD" = 0)) 
  trace <- 
    0:params$n_cycles %>% map_df(~({  
        apply(params$mP,3, simplify=FALSE, function(tp) {  
          (tr_ %*% (tp %^% .x)) %>% data.frame()
        })
    })) %>% 
    as.matrix() 
}

Specify the initial state occupancy.

Construct the Markov Trace

sim_cohort <- function(params) {
  tr_ <- t(c("Healthy" = params$n_cohort, "CVD" = 0)) 
  trace <- 
    0:params$n_cycles %>% map_df(~({  
        apply(params$mP,3, simplify=FALSE, function(tp) {  
          (tr_ %*% (tp %^% .x)) %>% data.frame()
        })
    })) %>% 
    as.matrix() 
}

Iterate over the cycles in the model.

Construct the Markov Trace

sim_cohort <- function(params) {
  tr_ <- t(c("Healthy" = params$n_cohort, "CVD" = 0)) 
  trace <- 
    0:params$n_cycles %>% map_df(~({  
        apply(params$mP,3, simplify=FALSE, function(tp) {  
          (tr_ %*% (tp %^% .x)) %>% data.frame()
        })
    })) %>% 
    as.matrix() 
}

Within a cycle, iterate over the strategies under consideration.

Construct the Markov Trace

The markov trace tells us the fraction of the population in each state in any cycle:

trace <- sim_cohort(params)
trace[1:5,]

     natural_history.Healthy natural_history.CVD
[1,]                 1.00000             0.00000
[2,]                 0.86071             0.13929
[3,]                 0.74082             0.25918
[4,]                 0.63763             0.36237
[5,]                 0.54881             0.45119

Checkpoint 1: CVD State Occupancy at Two Years

1. CVD State Occupancy at Two Years

Suppose we start off with a cohort of 1000 healthy individuals.
How many are in the CVD health state at year 2?

1. CVD State Occupancy at Two Years

Suppose we start off with a cohort of 1000 healthy individuals.
How many are in the CVD health state at year 2?

Pin this number to your brain for a few minutes!

(1000*trace[3,]) %>% t() %>% 
  data.frame() %>% 
  kable(col.names = c("Healthy","CVD")) %>% 
  kable_styling()

Healthy	CVD
740.82	259.18

Checkpoint 2: Life Expectancy

2. Life Expectancy

We have a model with 100 annual cycles and no death.
Let’s verify that life expectancy is 100.

2. Life Expectancy

We have a model with 100 annual cycles and no death.
Let’s verify that life expectancy is 100

1payoff_life_exp = c("Healthy" = 1, "CVD" = 1)
2life_exp_cycle = trace[-1,] %*% payoff_life_exp
3total_life_exp = sum(life_exp_cycle)

1: Life expectancy payoff is 1 for every state where someone remains alive.
2: Total life-years alive in each state is the markov trace multiplied by the payoff vector.
3: Total life expectancy is the sum of cycle values.

2. Life Expectancy

We have a model with 100 annual cycles and no death.
Let’s verify that life expectancy is 100

payoff_life_exp = c("Healthy" = 1, "CVD" = 1)
life_exp_cycle = trace[-1,] %*% payoff_life_exp
total_life_exp = sum(life_exp_cycle)
total_life_exp * params$cycle

[1] 100

Checkpoint 3: Changing Cycle Lengths

3. Daily Cycle Length

We will next switch to a daily cycle length and verify that we can recapitulate these numbers.

3. Daily Cycle Length

We will next switch to a daily cycle length and verify that we can recapitulate these numbers.

params_daily <- modifyList(params,list(
  cycle = 1/365,
  n_cycles = params$n_cycles * 365))

3. Daily Cycle Length

params_daily$mP <- 
  with(params_daily,{
    tp_H_CVD = 1 - exp(-r_H_CVD * cycle) 
    array(data = c(1 - tp_H_CVD, 0, 
                 tp_H_CVD,1),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
  })

apply(params_daily$mP,3,function(x) x, simplify=FALSE)

$natural_history
         to
from      Healthy        CVD
  Healthy 0.99959 0.00041087
  CVD     0.00000 1.00000000

3. Daily Cycle Length

# Markov Trace (Daily Cycle)
trace_daily <- 
  sim_cohort(params_daily)

# State Occupancy at 2 Years
1000 * trace_daily[365 * 2 + 1,]

trace_daily <- 
  sim_cohort(params_daily)
(1000 * trace_daily[365 * 2 + 1,]) %>% t() %>% 
data.frame() %>% 
  kable(col.names = c("Healthy","CVD")) %>% 
  kable_styling()

Healthy	CVD
740.82	259.18

3. Daily Cycle Length

# Life Expectancy
payoff_life_exp = c("Healthy" = 1, "CVD" = 1) 
life_exp_cycle_daily = trace_daily[-1,] %*% payoff_life_exp  
total_life_exp_daily = sum(life_exp_cycle_daily) 
total_life_exp_daily * params_daily$cycle

[1] 100

Summary

After 2 years, 259 individuals are in the CVD health state.
Because there is no death, after 100 cycles life expectancy is 100 years.
We verified that the above results hold exactly when we convert to a daily cycle length.

Part 2: Include Background Mortality

We next include a background mortality rate that is increased by a hazard ratio (10.0) among people with CVD.

We next include a background mortality rate that is increased by a hazard ratio (10.0) among people with CVD.

Parameterize the model

params_mort = 
  list(
    t_names = c("natural_history"),           # Strategy names. 
    n_treatments = 1,                         # Number of treatments
    s_names  = c("Healthy", "CVD", "Dead"),   # State names
    n_states = 3,                             # Number of states
    n_cohort = 1,                             # Cohort size
    n_cycles = 100,                           # Number of cycles in model.  
    cycle = 1,                                # Cycle length
    initial_age = 55,                         # Cohort starting age
    r_H_CVD = 0.15,                           # Rate of healthy -> CVD
    hr_CVD = 10,                              # Hazard Ratio: CVD Death
    r_H_D = 0.01                              # Rate of healthy -> dead
  )

Transition Probability Matrix

params_mort$mP <- 
  with(params_mort,{
    
    tp_H_CVD = 1 - exp(-r_H_CVD * cycle) 
    tp_H_D = 1 - exp(-r_H_D * cycle)
    tp_CVD_D <- 1 - exp(-(hr_CVD * r_H_D) * cycle)
    
    array(data = c(1 - tp_H_CVD - tp_H_D, 0, 0, 
                 tp_H_CVD,1-tp_CVD_D, 0,
                 tp_H_D, tp_CVD_D, 1),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
  })

Transition Probability Matrix

$natural_history
         to
from      Healthy     CVD      Dead
  Healthy 0.85076 0.13929 0.0099502
  CVD     0.00000 0.90484 0.0951626
  Dead    0.00000 0.00000 1.0000000

Construct the Markov Trace

sim_cohort <- function(params) {
  tr_ <- c(1,rep(0,params$n_states-1))
  trace <- 
    0:params$n_cycles %>% map_df(~({ 
        apply(params$mP,3, simplify=FALSE, function(tp) {  
          (tr_ %*% (tp %^% .x)) %>% data.frame()
        })
    })) %>% 
    as.matrix() 
}

trace_mort <- 
  sim_cohort(params_mort)

Construct the Markov Trace

trace_mort[1:5,]

     natural_history.Healthy natural_history.CVD natural_history.Dead
[1,]                 1.00000             0.00000            0.0000000
[2,]                 0.85076             0.13929            0.0099502
[3,]                 0.72379             0.24454            0.0316707
[4,]                 0.61577             0.32209            0.0621437
[5,]                 0.52387             0.37721            0.0989213

1. CVD State Occupancy at 2 Years

(1000*trace_mort[3,]) %>% t() %>% 
  data.frame() %>% 
  kable(digits = 1,col.names = c("Healthy","CVD","Dead")) %>% 
  kable_styling()

Healthy	CVD	Dead
723.8	244.5	31.7

2. Life Expectancy

payoff_life_exp_mort = c("Healthy" = 1, "CVD" = 1, "Dead" = 0) 
life_exp_cycle_mort = trace_mort[-1,] %*% payoff_life_exp_mort 
total_life_exp_mort = sum(life_exp_cycle_mort)
total_life_exp_mort

[1] 15.507

3. Daily Cycle Length

params_mort_daily <- modifyList(params_mort,list(
  cycle = 1/365,
  n_cycles = params$n_cycles * 365))

3. Daily Cycle Length

params_mort_daily$mP <- 
  with(params_mort_daily,{
    tp_H_CVD = 1 - exp(-r_H_CVD * cycle) 
    tp_H_D = 1 - exp(-r_H_D * cycle)
    tp_CVD_D <- 1 - exp(-(hr_CVD * r_H_D) * cycle)
    
    array(data = c(1 - tp_H_CVD - tp_H_D, 0, 0, 
                 tp_H_CVD,1-tp_CVD_D, 0,
                 tp_H_D, tp_CVD_D, 1),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
  })

3. Daily Cycle Length

apply(params_mort_daily$mP,3,function(x) x, simplify = FALSE)

$natural_history
         to
from      Healthy        CVD        Dead
  Healthy 0.99956 0.00041087 0.000027397
  CVD     0.00000 0.99972606 0.000273935
  Dead    0.00000 0.00000000 1.000000000

3. Daily Cycle Length

# Markov Trace (Daily Cycle)
trace_mort_daily <- sim_cohort(params_mort_daily)

Healthy	CVD	Dead
726.1	231.5	42.4

3. Daily Cycle Length

# Life Expectancy
payoff_mort_life_exp = c("Healthy" = 1, "CVD" = 1, "Death" = 0) 
life_exp_mort_daily = trace_mort_daily[-1,] %*% payoff_mort_life_exp  
total_life_exp_mort_daily = sum(life_exp_mort_daily) 
total_life_exp_mort_daily * params_mort_daily$cycle

[1] 15.624

Signs of Trouble

In the presence of competing events (mortality), we obtain different results for different cycle lengths!

Signs of Trouble

In the presence of competing events (mortality), we obtain different results for different cycle lengths!

cvd2_ann <- unname((1000*trace[3,2]))
cvd2_day <- unname((1000*trace_daily[365*2+1,2]))
cvd2_ann_mort <- 1000*trace_mort[3,2]
cvd2_day_mort <- 1000*trace_mort_daily[365*2+1,2]

le_ann <- total_life_exp * params$cycle
le_day <- total_life_exp_daily * params_daily$cycle
le_mort_ann <- total_life_exp_mort * params_mort$cycle
le_mort_day <- total_life_exp_mort_daily * params_mort_daily$cycle 

tibble(measure = c("1. CVD State Occupancy at Two Years", "2. Life Expectancy"), 
       annual = c(cvd2_ann, le_ann),
       daily = c(cvd2_day, le_day),
       annual_mort= c(cvd2_ann_mort,le_mort_ann),
       daily_mort = c(cvd2_day_mort, le_mort_day)) %>% 
  gt()  %>% 
  gt::data_color(
    method = "numeric",
    palette = "ggsci::red_material",
    domain = c(15.5,15.65),
    columns = c(annual_mort, daily_mort),
    rows = measure == "2. Life Expectancy"
  ) %>% 
  gt::data_color(
    method = "numeric",
    palette = "ggsci::red_material",
    domain = c(230,245),
    columns = c(annual_mort, daily_mort),
    rows = measure == "1. CVD State Occupancy at Two Years"
  ) %>% 
  tab_spanner(label = "Model With\nNo Mortality", columns = c(annual,daily)) %>% 
  tab_spanner(label = "Model With Mortality", columns = c(annual_mort, daily_mort)) %>% 
  cols_label(
    annual = "Yearly Cycle",
    daily = "Daily Cycle",
    annual_mort = "Yearly Cycle", 
    daily_mort = "Daily Cycle",
    measure = ""
  ) %>% 

  rm_stubhead() %>% 
  tab_options(
    table.font.size = pct(90)
    ) %>% 
   as_raw_html()

	Model With No Mortality		Model With Mortality
	Yearly Cycle	Daily Cycle	Yearly Cycle	Daily Cycle
1. CVD State Occupancy at Two Years	259.18	259.18	244.540	231.488
2. Life Expectancy	100.00	100.00	15.507	15.624

Signs of Trouble

In the presence of competing events (mortality), we obtain different results for different cycle lengths! Notice how the state occupancy in CVD differs in the right panel of the figure below.

df_res <- 
    data.frame(trace) %>% 
        mutate(cycle = (row_number()-1)) %>% 
        mutate(strategy = "Model With\nNo Mortality") %>% 
        mutate(cycle = cycle * params$cycle) %>% 
        mutate(cycle_length = paste0(params$cycle)) %>% 
    bind_rows(
        data.frame(trace_daily) %>% 
        mutate(cycle = (row_number()-1)) %>% 
        mutate(strategy = "Model With\nNo Mortality") %>% 
        mutate(cycle = cycle * params_daily$cycle) %>% 
        mutate(cycle_length = paste0(params_daily$cycle))
    ) %>% 
    bind_rows(
        data.frame(trace_mort) %>% 
        mutate(cycle = (row_number()-1)) %>% 
        mutate(strategy = "Model With Mortality") %>% 
        mutate(cycle = cycle * params_mort$cycle) %>% 
        mutate(cycle_length = paste0(params_mort$cycle))
    ) %>%  
    bind_rows(
        data.frame(trace_mort_daily) %>% 
        mutate(cycle = (row_number()-1)) %>% 
        mutate(strategy = "Model With Mortality") %>% 
        mutate(cycle = cycle * params_mort_daily$cycle) %>% 
        mutate(cycle_length = paste0(params_mort_daily$cycle))
    ) %>%  
    mutate(cycle_length = factor(round(as.numeric(paste0(cycle_length)),2), labels = c("Day","Year"))) %>% 
    gather(state,value,-strategy,-cycle,-cycle_length)  %>% 
    mutate(state = gsub("natural_history.","",state)) %>% 
   mutate(strategy = factor(strategy, levels = c("Model With\nNo Mortality","Model With Mortality")))

df_res %>% 
    ggplot(aes(x = cycle, y = value, colour = state ,lty=cycle_length)) + 
    geom_line() +
    facet_grid(~strategy) +
    hrbrthemes::theme_ipsum(base_family = "Arial") + 
    ggsci::scale_color_aaas(name="") + 
    theme(legend.position = "top") +
  # TK ADD BACK IN
    geom_dl(method = list("last.bumpup",hjust=1,vjust=-1),aes(label = state)) +
    geom_dl(method = list("first.bumpup",hjust=0,vjust=-1),aes(label = state)) + 
    guides(color = "none") + 
  scale_linetype_discrete(name="Cycle Length")

Warning: Removed 36602 rows containing missing values (`geom_line()`).

Warning: Removed 36602 rows containing missing values (`geom_dl()`).
Removed 36602 rows containing missing values (`geom_dl()`).

What Went Wrong?

To gain intuition for what is going on, let’s split out death for those with vs. without CVD:

What Went Wrong?

Parameterize the model

params_mort2 = 
  list(
    t_names = c("natural_history"),                       
    n_treatments = 1,                                     
    s_names  = c("Healthy", "CVD", "CVDDeath", "Dead"),   
    n_states = 4,                                         
    n_cohort = 1,                                         
    n_cycles = 100,                                        
    cycle = 1,                                            
    initial_age = 55,                                     
    r_H_CVD = 0.15,                                       
    hr_CVD = 10,                                          
    r_H_D = 0.01                                          
  )

params_mort2$mP <- 
  with(params_mort2,{
    
    tp_H_CVD = 1 - exp(-r_H_CVD * cycle) 
    tp_H_D = 1 - exp(-r_H_D * cycle)
    tp_CVD_D <- 1 - exp(-(hr_CVD * r_H_D) * cycle)
    
    array(data = c(1 - tp_H_CVD - tp_H_D, 0, 0,0,  
                 tp_H_CVD,1-tp_CVD_D, 0,0,
                 0,tp_CVD_D,1,0,
                 tp_H_D, 0,0, 1),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
  })

Spot the Problem (1yr Cycle)

apply(params_mort2$mP,3,function(x) x, simplify=FALSE)

$natural_history
          to
from       Healthy     CVD CVDDeath      Dead
  Healthy  0.85076 0.13929 0.000000 0.0099502
  CVD      0.00000 0.90484 0.095163 0.0000000
  CVDDeath 0.00000 0.00000 1.000000 0.0000000
  Dead     0.00000 0.00000 0.000000 1.0000000

Spot the Problem (1yr Cycle)

apply(params_mort2$mP,3,function(x) x, simplify=FALSE)

$natural_history
          to
from       Healthy     CVD CVDDeath      Dead
  Healthy  0.85076 0.13929 0.000000 0.0099502
  CVD      0.00000 0.90484 0.095163 0.0000000
  CVDDeath 0.00000 0.00000 1.000000 0.0000000
  Dead     0.00000 0.00000 0.000000 1.0000000

Recall that in this model, CVD has a very high death rate—meaning that many people who develop CVD will die very soon after.

Spot the Problem (1yr Cycle)

apply(params_mort2$mP,3,function(x) x, simplify=FALSE)

$natural_history
          to
from       Healthy     CVD CVDDeath      Dead
  Healthy  0.85076 0.13929 0.000000 0.0099502
  CVD      0.00000 0.90484 0.095163 0.0000000
  CVDDeath 0.00000 0.00000 1.000000 0.0000000
  Dead     0.00000 0.00000 0.000000 1.0000000

Recall that in this model, CVD has a very high death rate—meaning that many people who develop CVD will die very soon after.

With this knowledge in mind, can you spot something wrong with the above (annual) transition matrix?

Spot the Problem (10yr Cycle)

params_mort_10y = 
  list(
    t_names = c("natural_history"),           # Strategy names. 
    n_treatments = 1,                         # Number of treatments
    s_names  = c("Healthy", "CVD", "Dead"),   # State names
    n_states = 3,                             # Number of states
    n_cohort = 1,                             # Cohort size
    n_cycles = 10,                           # Number of cycles in model.  
    cycle = 10,                                # Cycle length
    initial_age = 55,                         # Cohort starting age
    r_H_CVD = 0.15,                           # Rate of healthy -> CVD
    hr_CVD = 10,                              # Hazard Ratio: CVD Death
    r_H_D = 0.01                              # Rate of healthy -> dead
  )

params_mort_10y$mP <- 
  with(params_mort_10y,{
    
    tp_H_CVD = 1 - exp(-r_H_CVD * cycle) 
    tp_H_D = 1 - exp(-r_H_D * cycle)
    tp_CVD_D <- 1 - exp(-(hr_CVD * r_H_D) * cycle)
    
    array(data = c(1 - tp_H_CVD - tp_H_D, 0, 0, 
                 tp_H_CVD,1-tp_CVD_D, 0,
                 tp_H_D, tp_CVD_D, 1),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
  })


params_mort2_10y = 
  list(
    t_names = c("natural_history"),                       # Strategy names. 
    n_treatments = 1,                                     # Number of treatments
    s_names  = c("Healthy", "CVD", "CVDDeath", "Dead"),   # State names
    n_states = 4,                                         # Number of states
    n_cohort = 1,                                         # Cohort size
    n_cycles = 10,                                       # Number of cycles in model.  
    cycle = 10,                                            # Cycle length
    initial_age = 55,                                     # Cohort starting age
    r_H_CVD = 0.15,                                       # Rate of healthy -> CVD
    hr_CVD = 10,                                          # Hazard Ratio: CVD Death
    r_H_D = 0.01                                          # Rate of healthy -> dead
  )

params_mort2_10y$mP <- 
  with(params_mort2_10y,{
    
    tp_H_CVD = 1 - exp(-r_H_CVD * cycle) 
    tp_H_D = 1 - exp(-r_H_D * cycle)
    tp_CVD_D <- 1 - exp(-(hr_CVD * r_H_D) * cycle)
    
    array(data = c(1 - tp_H_CVD - tp_H_D, 0, 0,0,  
                 tp_H_CVD,1-tp_CVD_D, 0,0,
                 0,tp_CVD_D,1,0,
                 tp_H_D, 0,0, 1),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
  })
apply(params_mort2_10y$mP,3,function(x) x, simplify = FALSE)

$natural_history
          to
from       Healthy     CVD CVDDeath     Dead
  Healthy  0.12797 0.77687  0.00000 0.095163
  CVD      0.00000 0.36788  0.63212 0.000000
  CVDDeath 0.00000 0.00000  1.00000 0.000000
  Dead     0.00000 0.00000  0.00000 1.000000

Recall that in this model, CVD has a very high death rate—meaning that many people who develop CVD will die very soon after.

With this knowledge in mind, can you spot something wrong with the above (annual) transition matrix?

Discretizing a Continuous Process

Our constructed model rules out compound transitions within a cycle.
Compound transition: Healthy CVD CVD Death in a single cycle.

apply(params_mort2_10y$mP,3,function(x) x, simplify = FALSE)

$natural_history
          to
from       Healthy     CVD CVDDeath     Dead
  Healthy  0.12797 0.77687  0.00000 0.095163
  CVD      0.00000 0.36788  0.63212 0.000000
  CVDDeath 0.00000 0.00000  1.00000 0.000000
  Dead     0.00000 0.00000  0.00000 1.000000

Discretizing a Continuous Process

To accurately capture the underlying continuous time process, there should be a transition probability defined for a seemingly impossible transition: Healthy CVD Death.

apply(params_mort2_10y$mP,3,function(x) x, simplify = FALSE)

$natural_history
          to
from       Healthy     CVD CVDDeath     Dead
  Healthy  0.12797 0.77687  0.00000 0.095163
  CVD      0.00000 0.36788  0.63212 0.000000
  CVDDeath 0.00000 0.00000  1.00000 0.000000
  Dead     0.00000 0.00000  0.00000 1.000000

Discretizing a Continuous Process

Why Don’t the Conversion Formulas Work?

When we use the standard rate-to-probability conversion formula (i.e., \(p_{H,CVD}= 1 - \exp(-r_{H,CVD} t)\) ) we obtain the marginal transition probability.
This marginal probability reflects the union of all possible transitions that start Healthy and then transition to and from the CVD state in the cycle.

Why Don’t the Conversion Formulas Work?

That is, \(p_{H,CVD}= 1 - \exp(-r_{H,CVD} t)\) captures the probability of following scenarios occurring in a cycle:

Individual transitions from healthy to CVD.
Individual transitions from healthy to CVD to CVD Death.
Individual dies from background causes before they would have developed CVD in the same cycle.

10 Year Cycle

apply(params_mort2_10y$mP,3,function(x) x, simplify = FALSE)

$natural_history
          to
from       Healthy     CVD CVDDeath     Dead
  Healthy  0.12797 0.77687  0.00000 0.095163
  CVD      0.00000 0.36788  0.63212 0.000000
  CVDDeath 0.00000 0.00000  1.00000 0.000000
  Dead     0.00000 0.00000  0.00000 1.000000

This matrix “captures” all transitions from Healthy CVD in a single transition probability (0.7769).

10 Year Cycle

apply(params_mort2_10y$mP,3,function(x) x, simplify = FALSE)

$natural_history
          to
from       Healthy     CVD CVDDeath     Dead
  Healthy  0.12797 0.77687  0.00000 0.095163
  CVD      0.00000 0.36788  0.63212 0.000000
  CVDDeath 0.00000 0.00000  1.00000 0.000000
  Dead     0.00000 0.00000  0.00000 1.000000

But we actually need to distribute out the marginal probability so that we’re left only with the probability that someone transitions to CVD and remains there at the end of the cycle.
As constructed, we have ruled out the possibility of compound transitions (Healthy CVD CVDDeath) in a single cycle.

10 Year Cycle

apply(params_mort_10y$mP,3,function(x) x, simplify = FALSE)

$natural_history
         to
from      Healthy     CVD     Dead
  Healthy 0.12797 0.77687 0.095163
  CVD     0.00000 0.36788 0.632121
  Dead    0.00000 0.00000 1.000000

This problem isn’t fixed in the version of the model where death is a single health state.
The problem persists, but is “hidden” from view!

10 Year Cycle

apply(params_mort2_10y$mP,3,function(x) x, simplify = FALSE)

$natural_history
          to
from       Healthy     CVD CVDDeath     Dead
  Healthy  0.12797 0.77687  0.00000 0.095163
  CVD      0.00000 0.36788  0.63212 0.000000
  CVDDeath 0.00000 0.00000  1.00000 0.000000
  Dead     0.00000 0.00000  0.00000 1.000000

This results in overcounting CVD state occupancy in some cases (people who develop CVD and die soon after) and undercounting CVD state occupancy in others (people who die of background causes before they would have developed CVD in the same cycle).

Discretizing a Continuous Process

Changing to shorter cycle does not get us out of this problem.

Quantified Absolute Error

Error after one cycle

truth <- Vectorize(function(r1, r2, t)
{
  # Embed properly
  x <- matrix(c(-r1-r2, r1, r2, rep(0, 6)), nrow=3, ncol=3, byrow=TRUE)
  (c(1, 0, 0) %*% expm(x*t))[,2:3]
})

misspecify <- Vectorize(function(r1, r2, t)
{
  # Common Advice
  # RL Fleurence, Christopher S. Hollenbeak. Rates and Probabilities in Economics, Pharmacoeconomics 2007; 25 (1) 3-6
  tp1 <- 1-exp(-r1*t)
  tp2 <- 1-exp(-r2*t)
 
  # THIS CONVERSION CAN RESULT IN IMPOSSIBLE VALUES!
  if(tp1 + tp2 > 1) return(NA)
  
  # Mispecify
  x <- matrix(c(1-tp1-tp2, tp1, tp2, 0, 1, 0, 0, 0, 1), nrow=3, ncol=3, byrow=TRUE)
  (c(1, 0, 0) %*% x)[,2:3]
})

par(mfrow=c(1,1))
abs.error <- Vectorize(function(r1, r2, t=1)
{
  misspecify(r1, r2, t) - truth(r1, r2, t)
})

x <- seq(0.01, 0.7, by = 0.01)
y <- x
z <- outer(x, y, FUN=function(x, y) {abs.error(x, y)[1,]})
contour(x,y,z, xlab="Rate 1 A -> B", ylab="Rate 2 A -> C",
        main="Misspecified Competing Rates",
        sub="Absolute Error Contour State B, t=1")

Quantified Relative Error

rel.error <- Vectorize(function(r1, r2, t=1)
{
    x <- truth(r1, r2, t)
    100*(misspecify(r1, r2, t) - x)/x
})

par(mfrow=c(1,2))
z <- outer(x, y, FUN=function(x, y) {rel.error(x, y)[1,]})
contour(x,y,z, xlab="Rate 1 A -> B", ylab="Rate 2 A -> C",
        main="",
        sub="Relative Error Contour (%) State B, t=1")

z <- outer(x, y, FUN=function(x, y) {rel.error(x, y)[2,]})
contour(x,y,z, xlab="Rate 1 A -> B", ylab="Rate 2 A -> C",
        main="",
        sub="Relative Error Contour (%) State C, t=1")
par(mfrow=c(1,1))
title(main="Misspecified Competing Rates")

For rates < 0.1 of the chosen time step the error is <5%.
This error will accumulate over each timestep.

Comparisons Across Strategies Bail You Out

Common CEA approaches compare strategies, i.e., \(f(\theta_A) - f(\theta_{B})\)
- Net monetary benefit.
- Net health benefit.
- Incremental cost effectiveness ratio.
With error it becomes:

Comparisons Across Strategies Bail You Out

Common CEA approaches compare strategies, i.e., \(f(\theta_A) - f(\theta_{B})\)
- Net monetary benefit.
- Net health benefit.
- Incremental cost effectiveness ratio.
With error it becomes:

\[f(\theta_A) + \epsilon_A - (f(\theta_{B}) + \epsilon_{B})\]

Specification Error

If \(\epsilon_A - \epsilon_{B} \sim 0\) the impact to the decision threshold is minimal.

These decision thresholds are semi-robust against this misspecification.
Restated: Single model runs are biased under misspecification, but comparison across strategies helps cancel this out.

Rate Misspecification Conclusions

In the presence of competing events (mortality), using standard conversion formulas will yield incorrect transition probability matrices.

Why Does This Matter?

Without properly embedding, you are no longer modeling a continuous time process.
So your model results will differ compared to a continuous time model (e.g., discrete event simulation)
Microsimulation suffers the same problem because the patient-level transition probability matrix is embedded in a time step and requires Markov Rate conversion.

Rate Misspecification Conclusions

Existing publication’s decision thresholds are generally okay if they misspecified rates.
Check the maximum rate of an event as a rough estimate of potential bias.
Rule of thumb, on the error plot use the estimated rate for full time scale of simulation, i.e. one step is entire simulation.

Rate Misspecification Error Decision Threshold

The amount of error for a decision threshold is the difference in contour. In this example: \(0.07 - 0.02 = 0.05\).

truth <- Vectorize(function(r1, r2, t)
{
  # Embed properly
  x <- matrix(c(-r1-r2, r1, r2, rep(0, 6)), nrow=3, ncol=3, byrow=TRUE)
  (c(1, 0, 0) %*% expm(x*t))[,2:3]
})

misspecify <- Vectorize(function(r1, r2, t)
{
  # Common Advice
  # RL Fleurence, Christopher S. Hollenbeak. Rates and Probabilities in Economics, Pharmacoeconomics 2007; 25 (1) 3-6
  tp1 <- 1-exp(-r1*t)
  tp2 <- 1-exp(-r2*t)
 
  # THIS CONVERSION CAN RESULT IN IMPOSSIBLE VALUES!
  if(tp1 + tp2 > 1) return(NA)
  
  # Mispecify
  x <- matrix(c(1-tp1-tp2, tp1, tp2, 0, 1, 0, 0, 0, 1), nrow=3, ncol=3, byrow=TRUE)
  (c(1, 0, 0) %*% x)[,2:3]
})

par(mfrow=c(1,1))
abs.error <- Vectorize(function(r1, r2, t=1)
{
  misspecify(r1, r2, t) - truth(r1, r2, t)
})

x <- seq(0.01, 0.7, by = 0.01)
y <- x
z <- outer(x, y, FUN=function(x, y) {abs.error(x, y)[1,]})
contour(x,y,z, xlab="Rate 1 A -> B", ylab="Rate 2 A -> C",
        main="Misspecified Competing Rates",
        sub="Absolute Error Contour State B, t=1")

points(x[8], y[63], col='red', lwd=4, cex=2)
points(x[63], y[37], col='red', lwd=4, cex=2)

lines(c(x[63], x[63]), c(y[37], y[9]), lty=2, lwd=2)

Additional Notes

Discrete Event Simulation automatically handles competing events.
Microsimulation misspecification error interacts with the discrete time step truncation and increases this error!
Published rates from clinical studies are usually adjusted for competing events (the purpose of survival studies).

So How Do We Obtain the Correct Transition Matrix?

Rather than use the conversion formulas, we need to properly “embed” the transition probability matrix.
The correct transition probability matrix includes a seemingly impossible transition: Healthy CVDDeath

Rather than use the conversion formulas, we need to properly “embed” the transition probability matrix.
The correct transition probability matrix includes a seemingly impossible transition: Healthy CVDDeath

params_mort2_10y_emb <- params_mort2_10y

params_mort2_10y_emb$mR <- 
  with(params_mort2_10y,{

    r_CVD_D <- hr_CVD * r_H_D
    
    R_ <- 
      array(data = c(0, 0, 0,0,  
                 r_H_CVD,0, 0,0,
                 0,r_CVD_D,0,0,
                 r_H_D, 0,0, 0),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
    R <- apply(R_,3, simplify=FALSE,function(x) {
      diag(x) = -rowSums(x)
      x * cycle
    })
  R    
  })
params_mort2_10y_emb$mP <- lapply(params_mort2_10y_emb$mR,expm)

round(params_mort2_10y_emb$mP[[1]],3)

         Healthy   CVD CVDDeath Dead
Healthy    0.202 0.415    0.333 0.05
CVD        0.000 0.368    0.632 0.00
CVDDeath   0.000 0.000    1.000 0.00
Dead       0.000 0.000    0.000 1.00

params_mort2_emb <- params_mort2

params_mort2_emb$mR <- 
  with(params_mort2_emb,{

    r_CVD_D <- hr_CVD * r_H_D
    
    R_ <- 
      array(data = c(0, 0, 0,0,  
                 r_H_CVD,0, 0,0,
                 0,r_CVD_D,0,0,
                 r_H_D, 0,0, 0),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
    R <- apply(R_,3, simplify=FALSE,function(x) {
      diag(x) = -rowSums(x)
      x * cycle
    })
  R    
  })
params_mort2_emb$mP <- lapply(params_mort2_emb$mR,expm)


sim_cohort_tmp <- function(params) {
  tr_ <- c(1,rep(0,params$n_states-1))
  trace <- 
    0:params$n_cycles %>% map_df(~({ 
        lapply(params$mP,function(tp) {  
          (tr_ %*% (tp %^% .x)) %>% data.frame()
        })
    })) %>% 
    as.matrix() 
  trace
}

trace_mort_emb <- 
  sim_cohort_tmp(params_mort2_emb)

trace_mort2 <- 
  sim_cohort(params_mort2)

# sum(trace_mort_emb[-1,] %*% c(1,1,0,0)) * params_mort2_emb$cycle
# sum(trace_mort2[-1,] %*% c(1,1,0,0)) * params_mort2$cycle

What Is Commonly Taught

Convert each rate to a probability for the time step using \(p = 1-\exp(-rt)\).
Place the probabilities within a transition matrix \(\mathbf{P}\).
Diagonal elements of \(\mathbf{P}\): 1 - (row sum of nondiagonal elements).

What We’ll Teach You

Place each rate in a rate matrix \(\mathbf{R}\).
Diagonal elements of \(\mathbf{R}\): -1 * (row sum of nondiagonal elements).
Embed the transition probability matrix using matrix exponentiation \(\mathbf{P} = \exp{(\mathbf{R}t)}\).

What We’ll Teach You

Place each rate in a rate matrix \(\mathbf{R}\).
Diagonal elements of \(\mathbf{R}\): -1 * (row sum of nondiagonal elements).
Embed the transition probability matrix using matrix exponentiation \(\mathbf{P} = \exp{(\mathbf{R}t})\).

This process will ensure that the matrix captures “jumpover” states for compound transitions!

What We’ll Teach You

Place each rate in a rate matrix \(\mathbf{R}\).
Diagonal elements of \(\mathbf{R}\): -1 * (row sum of nondiagonal elements).
Embed the transition probability matrix using matrix exponentiation \(\mathbf{P} = \exp{(\mathbf{R}t})\).

This process will ensure that the matrix captures “jumpover” states for compound transitions!

It also ensures that changing cycle lengths will not change your answers!

Steps to Embed a Transition Matrix

Parameterize the model

params_mort_embedded = 
  list(
    t_names = c("natural_history"),                       
    n_treatments = 1,                                     
    s_names  = c("Healthy", "CVD", "CVDDeath", "Dead"),   
    n_states = 4,                                         
    n_cohort = 1,                                         
    n_cycles = 100,                                        
    cycle = 1,                                            
    initial_age = 55,                                     
    r_H_CVD = 0.15,                                       
    hr_CVD = 10,                                          
    r_H_D = 0.01                                          
  )

Note: this is the same parameter list as earlier!

Transition Rate Matrix

params_mort_embedded$mR <- 
  with(params_mort2_emb,{

    r_CVD_D <- hr_CVD * r_H_D
    
    R_ <- 
      array(data = c(0, 0, 0,0,  
                 r_H_CVD,0, 0,0,
                 0,r_CVD_D,0,0,
                 r_H_D, 0,0, 0),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
    R <- apply(R_,3, simplify=FALSE,function(x) {
      diag(x) = -rowSums(x)
      x * cycle
    })
  R    
  })

Transition Rate Matrix

params_mort_embedded$mR <- 
  with(params_mort2_emb,{

    r_CVD_D <- hr_CVD * r_H_D
    
    R_ <- 
      array(data = c(0, 0, 0,0,  
                 r_H_CVD,0, 0,0,
                 0,r_CVD_D,0,0,
                 r_H_D, 0,0, 0),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
    R <- apply(R_,3, simplify=FALSE,function(x) {
      diag(x) = -rowSums(x)
      x * cycle
    })
  R    
  })

Calculate the CVD Dead rate

Transition Rate Matrix

params_mort_embedded$mR <- 
  with(params_mort2_emb,{

    r_CVD_D <- hr_CVD * r_H_D
    
    R_ <- 
      array(data = c(0, 0, 0,0,  
                 r_H_CVD,0, 0,0,
                 0,r_CVD_D,0,0,
                 r_H_D, 0,0, 0),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
    R <- apply(R_,3, simplify=FALSE,function(x) {
      diag(x) = -rowSums(x)
      x * cycle
    })
  R    
  })

Place the rates within the matrix. We’ll get the diagonal elements in a second to ensure things balance out.

Transition Rate Matrix

params_mort_embedded$mR <- 
  with(params_mort2_emb,{
    r_CVD_D <- hr_CVD * r_H_D
    R_ <- 
      array(data = c(0, 0, 0,0,  
                 r_H_CVD,0, 0,0,
                 0,r_CVD_D,0,0,
                 r_H_D, 0,0, 0),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
    R <- apply(R_,3, simplify=FALSE,function(x) {
      diag(x) = -rowSums(x)
      x * cycle 
    })
  R    
  })

Diagonal elements are the negative row sum (i.e., each row of \(\mathbf{R}\) sums to 0)

Transition Rate Matrix

params_mort_embedded$mR <- 
  with(params_mort2_emb,{
    r_CVD_D <- hr_CVD * r_H_D
    R_ <- 
      array(data = c(0, 0, 0,0,  
                 r_H_CVD,0, 0,0,
                 0,r_CVD_D,0,0,
                 r_H_D, 0,0, 0),
          dim = c(n_states, n_states, n_treatments),
                  dimnames = list(from = s_names,
                                  to = s_names,
                                  t_names))
    R <- apply(R_,3, simplify=FALSE,function(x) {
      diag(x) = -rowSums(x) 
      x * cycle
    })
  R    
  })

Convert \(\mathbf{R}\) to the model cycle length (e.g., 1 year, 1 day, 1 month, etc.).

Transition Rate Matrix

params_mort_embedded$mR

$natural_history
          to
from       Healthy   CVD CVDDeath Dead
  Healthy    -0.16  0.15      0.0 0.01
  CVD         0.00 -0.10      0.1 0.00
  CVDDeath    0.00  0.00      0.0 0.00
  Dead        0.00  0.00      0.0 0.00

Transition Probability Matrix

expm(params_mort_embedded$mR[["natural_history"]])

         Healthy     CVD  CVDDeath     Dead
Healthy  0.85214 0.13173 0.0068811 0.009241
CVD      0.00000 0.90484 0.0951626 0.000000
CVDDeath 0.00000 0.00000 1.0000000 0.000000
Dead     0.00000 0.00000 0.0000000 1.000000

Notice how the embedded matrix now has a transition probability from Healthy CVD Death!

Using Conversion Formulas:

round(apply(params_mort2$mP,3,function(x) x, simplify = FALSE)[[1]],3)

          to
from       Healthy   CVD CVDDeath Dead
  Healthy    0.851 0.139    0.000 0.01
  CVD        0.000 0.905    0.095 0.00
  CVDDeath   0.000 0.000    1.000 0.00
  Dead       0.000 0.000    0.000 1.00

Using Rate Matrix Eponentiation:

round(expm(params_mort_embedded$mR[["natural_history"]]),3)

         Healthy   CVD CVDDeath  Dead
Healthy    0.852 0.132    0.007 0.009
CVD        0.000 0.905    0.095 0.000
CVDDeath   0.000 0.000    1.000 0.000
Dead       0.000 0.000    0.000 1.000

Now it’s your turn!

Case Study

Learning Objectives

Construct a natural history model for CVD.
Background mortality based on cause-deleted life table data.
CVD mortality based on cause-specific mortality.
Markov model is properly embedded.
Markov model results will closely match cause-deleted and cause-specific mortality in the US population.